home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / plot_field.pro < prev    next >
Text File  |  1997-07-08  |  3KB  |  151 lines
  1. ; $Id: plot_field.pro,v 1.5 1996/12/17 21:12:53 ali Exp $
  2.  
  3. function mybi,a,x,y    ;Bilinear interpolation
  4. sizea=size(a)
  5. nx=sizea[1]
  6. i=long(x)+nx*long(y)
  7. q=y-long(y)
  8. p=x-long(x)
  9. aint=(1.0-p)*(1.0-q)*a[i]+p*(1.0-q)*a[i+nx]+q*(1.0-p)*a[i+1]+p*q*a[i+nx+1]
  10. return,aint
  11. end
  12.  
  13.  
  14.  
  15. PRO ARRHEAD,X        ;Add an arrowhead to a vector
  16. THETA=ATAN(1.0)/4.0
  17. TANT = TAN(THETA)
  18. NP=3.0
  19. SCAL=6.
  20.  
  21. SX=SIZE(X)
  22. N=SX[2]
  23.  
  24.  
  25. BIGL=SQRT((X[*,N-4,0]-X[*,N-5,0])^2+(X[*,N-4,1]-X[*,N-5,1])^2)
  26. wbigl=where(BIGL ne 0.0)
  27. wnbigl=where(bigl eq 0.0)
  28. LL  = SCAL*TANT*BIGL[wbigl]/NP
  29.  
  30. DX = LL*(X[wbigl,N-4,1]-X[wbigl,N-5,1])/BIGL[wbigl]
  31. DY = LL*(X[wbigl,N-4,0]-X[wbigl,N-5,0])/BIGL[wbigl]
  32.  
  33. XM = X[wbigl,N-4,0]-(SCAL-1)*(X[wbigl,N-4,0]-X[wbigl,N-5,0])/NP
  34. YM = X[wbigl,N-4,1]-(SCAL-1)*(X[wbigl,N-4,1]-X[wbigl,N-5,1])/NP
  35.  
  36. X[wbigl,N-3,0] = XM-DX
  37. X[wbigl,N-2,0] = X[wbigl,N-4,0]
  38. X[wbigl,N-1,0] = XM+DX
  39.  
  40. X[wbigl,N-3,1] = YM+DY
  41. X[wbigl,N-2,1] = X[wbigl,N-4,1]
  42. X[wbigl,N-1,1] = YM-DY
  43.  
  44. if n_elements(wnbigl) gt 1 then begin
  45. X[wnbigl,N-3,0] = x[wnbigl,n-4,0]
  46. X[wnbigl,N-2,0] = X[wnbigl,n-4,0]
  47. X[wnbigl,N-1,0] = X[wnbigl,n-4,0]
  48.  
  49. X[wnbigl,N-3,1] = X[wnbigl,N-4,1]
  50. X[wnbigl,N-2,1] = X[wnbigl,N-4,1]
  51. X[wnbigl,N-1,1] = X[wnbigl,N-4,1]
  52. end
  53. return
  54. END
  55.  
  56.  
  57.  
  58. function arrows,u,v,n,length,xmax,lmax
  59. su=size(u)
  60. nx=su[1]
  61. ny=su[2]
  62. if n_params(0) lt 6 then lmax=sqrt(max(u^2+v^2))
  63. if n_params(0) lt 5 then xmax=1.0
  64. lth=0.1*length/lmax
  65. xt=randomu(seed,n)
  66. yt=randomu(seed,n)
  67. x=fltarr(n,13,2)
  68. x[0,0,0]=xt[*]
  69. x[0,0,1]=yt[*]
  70. for i=1,9 do begin
  71.  xt[0]=(nx-1)*x[*,i-1,0]
  72.  yt[0]=(ny-1)*x[*,i-1,1]
  73.  ut=mybi(u,xt,yt)
  74.  vt=mybi(v,xt,yt)
  75.  x[0,i,0]=x[*,i-1,0]+ut*lth
  76.  x[0,i,1]=x[*,i-1,1]+vt*lth
  77. end
  78. ARRHEAD,X
  79. return,x<1.0>0.0
  80. end
  81.  
  82.  
  83.  
  84.  
  85. PRO Plot_field,U,V,N=N,LENGTH=length,ASPECT=aspect, title= title
  86. ;+
  87. ; NAME:
  88. ;    PLOT_FIELD
  89. ;
  90. ; PURPOSE:
  91. ;    This procedure plots a 2-dimensional field.
  92. ;
  93. ; CATEGORY:
  94. ;    Plotting, two-dimensional.
  95. ;
  96. ; CALLING SEQUENCE:
  97. ;    PLOT_FIELD, U, V
  98. ;
  99. ; INPUTS:
  100. ;    U:    The 2-dimensional array giving the field vector at each
  101. ;        point in the U[X] direction.
  102. ;
  103. ;    V:    The 2-dimensional array giving the field vector at each
  104. ;        point in the V[Y] direction.
  105. ;
  106. ; KEYWORD PARAMETERS:
  107. ;    N:    The number of arrows to draw. The default is 200.
  108. ;
  109. ;    LENGTH:    The length of the longest field vector expressed as a fraction
  110. ;        of the plotting area. The default is 0.1.
  111. ;
  112. ;    ASPECT:    The aspect ratio of the plot (i.e., the ratio of the X size 
  113. ;        to Y size). The default is 1.0.
  114. ;
  115. ;    TITLE:    The title of plot. The default is "Velocity Field".
  116. ;
  117. ; OUTPUTS:
  118. ;    No explicit outputs.
  119. ;
  120. ; SIDE EFFECTS:
  121. ;    A new plot is drawn to the current output device.
  122. ;
  123. ; PROCEDURE:
  124. ;    N random points are picked, and from each point a path is traced
  125. ;    along the field. The length of the path is proportional to "LENGTH" 
  126. ;    and the field vector magnitude.
  127. ;
  128. ; EXAMPLE:
  129. ;    X = FINDGEN(20, 20)        ; Create array X
  130. ;    Y = FINDGEN(20, 20)*3        ; Create array Y
  131. ;    PLOT_FIELD, X, Y        ; Create plot
  132. ;    
  133. ; MODIFICATION HISTORY:
  134. ;    Jan., 1988    Neal Hurlburt, University of Colorado.
  135. ;-
  136.  
  137.  
  138. if n_elements(n) le 0 then n=200
  139. if n_elements(length) le 0 then length=.1
  140. if n_elements(aspect) le 0 then aspect=1.0
  141. if n_elements(title) le 0 then title = 'Velocity Field'
  142.  
  143. X=ARROWS(U,V,N,LENGTH)
  144. IF ASPECT GT 1. THEN position=[0.20,(0.5-0.30/ASPECT),0.90,(0.5+0.40/ASPECT)]$
  145.  else position=[(0.5-0.30*ASPECT),0.20,(0.5+0.40*ASPECT),0.90]
  146. plot,[0,1,1,0,0],[0,0,1,1,0],title=title,/nodata, pos = position
  147. FOR I=0,N-1 DO OPLOT,X[I,*,0],X[I,*,1]
  148. RETURN
  149. end
  150.  
  151.